home *** CD-ROM | disk | FTP | other *** search
/ Mac Magazin/MacEasy 21 / Mac Magazin and MacEasy Magazine CD - Issue 21.iso / Wissenschaft & Technik / yorick12vr1-nofpu folder / include / slice3.i < prev    next >
Text File  |  1996-02-13  |  50KB  |  1,473 lines
  1. /*
  2.    SLICE3.I
  3.    find 2D slices of a 3D hexahedral mesh
  4.  
  5.    $Id$
  6.  */
  7. /*    Copyright (c) 1996.  The Regents of the University of California.
  8.                     All rights reserved.  */
  9.  
  10. /*
  11.  * Caveats:
  12.  * (A) Performance is reasonably good, but may still be a factor of
  13.  *     several slower than what could be achieved in compiled code.
  14.  * (B) Only a simple in-memory mesh model is implemented here.
  15.  *     However, hooks are supplied for more interesting possibilities,
  16.  *     such as a large binary file resident mesh data base.
  17.  * (C) There is a conceptual difficulty with _walk3 for the case
  18.  *     of a quad face all four of whose edges are cut by the slicing
  19.  *     plane.  This can only happen when two opposite corners are
  20.  *     above and the other two below the slicing plane.  There are
  21.  *     three possible ways to connect the four intersection points in
  22.  *     two pairs: (1) // (2) \\ and (3) X  There is a severe problem
  23.  *     with (1) and (2) in that a consistent decision must be made
  24.  *     when connecting the points on the two cells which share the
  25.  *     face - that is, each face must carry information on which way
  26.  *     it is triangulated.  For a regular 3D mesh, it is relatively
  27.  *     easy to come up with a consistent scheme for triangulating faces,
  28.  *     but for a general unstructured mesh, each face itself must carry
  29.  *     this information.  This presents a huge challenge for data flow,
  30.  *     which I don't believe is worthwhile.  Because the X choice is
  31.  *     unique, and I don't see why we shouldn't use it here.
  32.  *     For contouring routines, we reject the X choice on esthetic
  33.  *     grounds, and perhaps that will prove to be the case here as
  34.  *     well - but I believe we should try the simple way out first.
  35.  *     In this case, we are going to be filling these polygons with
  36.  *     a color representing a function value in the cell.  Since the
  37.  *     adjacent cells should have nearly the same values, the X-traced
  38.  *     polygons will have nearly the same color, and I doubt there will
  39.  *     be an esthetic problem.  Anyway, the slice3 implemented
  40.  *     below produces the unique X (bowtied) polygons, rather than
  41.  *     attempting to choose between // or \\ (non-bowtied) alternatives.
  42.  *     Besides, in the case of contours, the trivial alternating
  43.  *     triangulation scheme is just as bad esthetically as every
  44.  *     zone triangulated the same way!
  45.  */
  46.  
  47. /* ------------------------------------------------------------------------ */
  48.  
  49. func plane3(normal, point)
  50. /* DOCUMENT plane3(normal, point)
  51.          or plane3([nx,ny,nz], [px,py,pz])
  52.  
  53.      returns [nx,ny,nz,pp] for the specified plane.
  54.  
  55.    SEE ALSO: slice3, mesh3
  56.  */
  57. {
  58.   /* the normal doesn't really need to be normalized, but this
  59.    * has the desirable side effect of blowing up if normal==0 */
  60.   normal/= abs(normal(1),normal(2),normal(3));
  61.   return grow(normal,normal(+)*point(+));
  62. }
  63.  
  64. /* ------------------------------------------------------------------------ */
  65.  
  66. func mesh3(x,y,z,..)
  67. /* DOCUMENT mesh3(x,y,z)
  68.          or mesh3(x,y,z, f1,f2,...)
  69.      or mesh3(xyz, f1,f2,...)
  70.      or mesh3(nxnynz, dxdydz, x0y0z0, f1,f2,...)
  71.  
  72.      make mesh3 argument for slice3, xyz3, getv3, etc., functions.
  73.      X, Y, and Z are each 3D coordinate arrays.  The optional F1, F2,
  74.      etc. are 3D arrays of function values (e.g. density, temperature)
  75.      which have one less value along each dimension than the coordinate
  76.      arrays.  The "index" of each zone in the returned mesh3 is
  77.      the index in these cell-centered Fi arrays, so every index from
  78.      one through the total number of cells indicates one real cell.
  79.      The Fi arrays can also have the same dimensions as X, Y, or Z
  80.      in order to represent point-centered quantities.
  81.  
  82.      If X has four dimensions and the length of the first is 3, then
  83.      it is interpreted as XYZ (which is the quantity actually stored
  84.      in the returned cell list).
  85.  
  86.      If X is a vector of 3 integers, it is interpreted as [nx,ny,nz]
  87.      of a uniform 3D mesh, and the second and third arguments are
  88.      [dx,dy,dz] and [x0,y0,z0] respectively.  (DXDYDZ represent the
  89.      size of the entire mesh, not the size of one cell, and NXNYNZ are
  90.      the number of cells, not the number of points.)
  91.  
  92.    SEE ALSO: slice3, xyz3, getv3, getc3
  93.  */
  94. {
  95.   /* other sorts of meshes are possible -- a mesh which lives
  96.    * in a binary file is an obvious example -- which would need
  97.    * different workers for xyz3, getv3, getc3, and iterator3
  98.    * iterator3_rect may be more general than the other three;
  99.    * as long as the cell dimensions are the car of the list
  100.    * which is the 2nd car of m3, it will work */
  101.   virtuals= _lst(xyz3_rect, getv3_rect, getc3_rect, iterator3_rect);
  102.  
  103.   dims= dimsof(x);
  104.   if (dims(1)==4 && dims(2)==3 && min(dims)>=2) {
  105.     xyz= &double(x);
  106.     n= 2*(!is_void(y))+(!is_void(z));
  107.     dims= grow([3], dims(3:5));
  108.   } else if (dims(1)==1 && dims(2)==3 && structof(0+x)==long && min(x)>0 &&
  109.          dimsof(y)(1)==1 && dimsof(z)(1)==1 &&
  110.          numberof(y)==3 && numberof(z)==3) {
  111.     xyz= &double([y, z]);
  112.     dims= grow([3], 1+x);
  113.     n= 0;
  114.     _car, virtuals, 1, xyz3_unif;
  115.   } else {
  116.     if (dims(1)!=3 || min(dims)<2 ||
  117.     dimsof(y)(1)!=3 || anyof(dimsof(y)!=dims) ||
  118.     dimsof(z)(1)!=3 || anyof(dimsof(z)!=dims))
  119.       error, "X,Y,Z are not viable 3D coordinate mesh arrays";
  120.     xyz= array(0.0, 3, dimsof(x));
  121.     xyz(1,..)= x;
  122.     xyz(2,..)= y;
  123.     xyz(3,..)= z;
  124.     xyz= &xyz;
  125.     n= 0;
  126.   }
  127.  
  128.   dim_cell= dims;
  129.   dim_cell(2:4)-= 1;
  130.   m3= _lst(virtuals, _lst(dim_cell(2:4), *xyz));
  131.   last= _cdr(m3);
  132.  
  133.   while (n || more_args()) {
  134.     if (n) {
  135.       if (n&2) {
  136.     f= y;
  137.     n&= 1;
  138.       } else {
  139.     f= z;
  140.     n= 0;
  141.       }
  142.     } else {
  143.       f= next_arg();
  144.     }
  145.     if (dimsof(f)(1)!=3 ||
  146.     (anyof(dimsof(f)!=dims) && anyof(dimsof(f)!=dim_cell))) {
  147.       error, "F"+pr1(i)+" is not a viable 3D cell value";
  148.     }
  149.     _cdr, last, 1, _lst(f);
  150.     last= _cdr(last);
  151.   }
  152.  
  153.   return m3;
  154. }
  155.  
  156. /* ------------------------------------------------------------------------ */
  157.  
  158. /* Ways that a list of polygons can be extracted:
  159.  * Basic idea:
  160.  *   (1) At each *vertex* of the cell list, a function value is defined.
  161.  *       This is the "slicing function", perhaps the equation of a plane,
  162.  *       perhaps some other vertex-centered function.
  163.  *   (2) The slice3 routine returns a list of cells for which the
  164.  *       function value changes sign -- that is, cells for which some
  165.  *       vertices have positive function values, and some negative.
  166.  *       The function values and vertex coordinates are also returned.
  167.  *   (3) The slice3 routine computes the points along the *edges*
  168.  *       of each cell where the function value is zero (assuming linear
  169.  *       variation along each edge).  These points will be vertices of
  170.  *       the polygons.  The routine also sorts the vertices into cyclic
  171.  *       order.
  172.  *   (4) A "color function" can be used to assign a "color" or other
  173.  *       value to each polygon.  If this function depends only on the
  174.  *       coordinates of the polygon vertices (e.g.- 3D lighting), then
  175.  *       the calculation can be done elsewhere.  There are two other
  176.  *       possibilities:  The color function might be a cell-centered
  177.  *       quantity, or a vertex-centered quantity (like the slicing
  178.  *       function) on the mesh.  In these cases, slice3 already
  179.  *       has done much of the work, since it "knows" cell indices,
  180.  *       edge interpolation coefficients, and the like.
  181.  *
  182.  * There are two particularly important cases:
  183.  * (1) Slicing function is a plane, coloring function is either a
  184.  *     vertex or cell centered mesh function.  Coloring function
  185.  *     might also be a *function* of one or more of the predefined
  186.  *     mesh functions.  If you're eventually going to sweep the whole
  187.  *     mesh, you want to precalculate it, otherwise on-the-fly might
  188.  *     be better.
  189.  * (2) Slicing function is a vertex centered mesh function,
  190.  *     coloring function is 3D shading (deferred).
  191.  *
  192.  * fslice(m3, vertex_list)
  193.  * vertex_list_iterator(m3, vertex_list, mesh3)
  194.  * fcolor(m3, vertex_list, fslice_1, fslice_2)
  195.  *   the coloring function may need the value of fslice at the vertices
  196.  *   in order to compute the color values by interpolation
  197.  * two "edge functions": one to detect edges where sign of fslice changes,
  198.  *   second to interpolate for fcolor
  199.  *
  200.  * slice3(m3, fslice, &nverts, &xyzverts, <fcolor>)
  201.  */
  202.  
  203. func slice3(m3, fslice, &nverts, &xyzverts, fcolor, nointerp, value=)
  204. /* DOCUMENT slice3, m3, fslice, nverts, xyzverts
  205.          or color_values= slice3(m3, fslice, nverts, xyzverts, fcolor)
  206.          or color_values= slice3(m3, fslice, nverts, xyzverts, fcolor, 1)
  207.  
  208.      slice the 3D mesh M3 using the slicing function FSLICE, returning
  209.      the lists NVERTS and XYZVERTS.  NVERTS is the number of vertices
  210.      in each polygon of the slice, and XYZVERTS is the 3-by-sum(NVERTS)
  211.      list of polygon vertices.  If the FCOLOR argument is present, the
  212.      values of that coloring function on the polygons are returned as
  213.      the value of the slice3 function (numberof(color_values) ==
  214.      numberof(NVERTS) == number of polygons).
  215.  
  216.      If the slice function FSLICE is a function, it should be of the
  217.      form:
  218.         func fslice(m3, chunk)
  219.      returning a list of function values on the specified chunk of the
  220.      mesh m3.  The format of chunk depends on the type of m3 mesh, so
  221.      you should use only the other mesh functions xyz3 and getv3 which
  222.      take m3 and chunk as arguments.  The return value of fslice should
  223.      have the same dimensions as the return value of getv3; the return
  224.      value of xyz3 has an additional first dimension of length 3.
  225.  
  226.      If FSLICE is a list of 4 numbers, it is taken as a slicing plane
  227.      with the equation FSLICE(+:1:3)*xyz(+)-FSLICE(4), as returned by
  228.      plane3.
  229.  
  230.      If FSLICE is a single integer, the slice will be an isosurface for
  231.      the FSLICEth variable associated with the mesh M3.  In this case,
  232.      the keyword value= must also be present, representing the value
  233.      of that variable on the isosurface.
  234.  
  235.      If FCOLOR is nil, slice3 returns nil.  If you want to color the
  236.      polygons in a manner that depends only on their vertex coordinates
  237.      (e.g.- by a 3D shading calculation), use this mode.
  238.  
  239.      If FCOLOR is a function, it should be of the form:
  240.         func fcolor(m3, cells, l, u, fsl, fsu, ihist)
  241.      returning a list of function values on the specified cells of the
  242.      mesh m3.  The cells argument will be the list of cell indices in
  243.      m3 at which values are to be returned.  l, u, fsl, fsu, and ihist
  244.      are interpolation coefficients which can be used to interpolate
  245.      from vertex centered values to the required cell centered values,
  246.      ignoring the cells argument.  See getc3 source code.
  247.      The return values should always have dimsof(cells).
  248.  
  249.      If FCOLOR is a single integer, the slice will be an isosurface for
  250.      the FCOLORth variable associated with the mesh M3.
  251.  
  252.      If the optional argument after FCOLOR is non-nil and non-zero,
  253.      then the FCOLOR function is called with only two arguments:
  254.         func fcolor(m3, cells)
  255.  
  256.    SEE ALSO: mesh3, plane3, xyz3, getv3, getc3, slice2, plfp
  257.  */
  258. {
  259.   nverts= xyzverts= [];
  260.  
  261.   /* handle the different possibilities for fslice */
  262.   if (!is_func(fslice)) {
  263.     if (is_void(value) &&
  264.     dimsof(fslice)(1)==1 && numberof(fslice)==4) {
  265.       normal= fslice(1:3);
  266.       projection= fslice(4);
  267.       fslice= _plane_slicer;
  268.     } else if (dimsof(fslice)(1)==0 && structof(0+fslice)==long) {
  269.       if (is_void(value))
  270.     error, "value= keyword required when FSLICE is mesh variable";
  271.       iso_index= fslice;
  272.       fslice= _isosurface_slicer;
  273.     } else {
  274.       error, "illegal form of FSLICE argument, try help,slice3";
  275.     }
  276.   }
  277.  
  278.   /* will need cell list if fcolor function to be computed */
  279.   need_clist= !is_void(fcolor);
  280.  
  281.   /* test the different possibilities for fcolor */
  282.   if (need_clist && !is_func(fcolor)) {
  283.     if (dimsof(fcolor)(1)!=0 || structof(0+fcolor)!=long) {
  284.       error, "illegal form of FCOLOR argument, try help,slice3";
  285.     }
  286.   }
  287.  
  288.   /* chunk up the m3 mesh and evaluate the slicing function to
  289.    * find those cells cut by fslice==0
  290.    * chunking avoids potentially disasterously large temporaries
  291.    */
  292.   _xyz3_save= 1;      /* flag for xyz3 */
  293.   got_xyz= 0;
  294.   ntotal= nchunk= 0;
  295.   results= [];
  296.   for (chunk=iterator3(m3) ;
  297.        !is_void(chunk) ;
  298.        chunk=iterator3(m3,chunk)) {
  299.  
  300.     /* get the values of the slicing function at the vertices of
  301.      * this chunk
  302.      */
  303.     _xyz3= [];        /* if fslice calls xyz3, xyz saved here */
  304.     fs= fslice(m3, chunk);
  305.  
  306.     /* will need cell list if fslice did not compute xyz */
  307.     got_xyz= !is_void(_xyz3);
  308.     need_clist|= !got_xyz;
  309.  
  310.     /* If the m3 mesh is totally unstructured, the chunk should be
  311.      * arranged so that fslice returns an 2-by-2-by-2-by-ncells array
  312.      * of vertex values of the slicing function.
  313.      * On the other hand, if the mesh vertices are arranged in a
  314.      * rectangular grid (or a few patches of rectangular grids), the
  315.      * chunk should be the far less redundant rectangular patch.
  316.      */
  317.     dims= dimsof(fs);
  318.     irregular= numberof(dims)>4;
  319.     if (irregular) {
  320.       /* fs is a 2-by-2-by-2-by-ncells array
  321.        * here is the fastest way to generate the required cell list
  322.        */
  323.       clist= where((fs<0.0)(sum:1:8,1,1,) & 7);
  324.  
  325.     } else {
  326.       /* fs is an ni-by-nj-by-nk array
  327.        * result of the zcen is 0, 1/8, 2/8, ..., 7/8, or 1
  328.        */
  329.       clist= double(fs<0.0)(zcen,zcen,zcen);
  330.       clist= where(clist>.1 & clist<.9);
  331.       /* alternative possibilities which run about equally quickly are:
  332.        *    clist!=floor(clist)
  333.        *    clist%1.0
  334.        * these both rely on the fact that 0.5*(1.0+1.0)==1.0 exactly
  335.        * they also both call slow libm functions (floor and amod)
  336.        */
  337.     }
  338.  
  339.     if (numberof(clist)) {
  340.       ntotal+= numberof(clist);
  341.       /* we need to save:
  342.        * (1) the absolute cell indices of the cells in clist
  343.        * (2) the corresponding 2-by-2-by-2-by-ncells list of fslice
  344.        *     values at the vertices of these cells
  345.        */
  346.       if (irregular) {
  347.     fs= fs(,,,clist);
  348.     if (got_xyz) _xyz3= _xyz3(,,,,clist);
  349.       } else {
  350.     list= to_corners3(clist, dims(2), dims(3));
  351.     fs= fs(list);
  352.     if (got_xyz) _xyz3= _xyz3(,list);
  353.       }
  354.       /* here, the iterator converts to absolute cell indices without
  355.        * incrementing the chunk */
  356.       if (need_clist) clist= iterator3(m3, chunk, clist);
  357.       else clist= [];
  358.       if (!(nchunk&127)) grow, results, array(pointer, 3, nchunk+128);
  359.       ++nchunk;
  360.       results(1,nchunk)= &clist;
  361.       results(2,nchunk)= &fs;
  362.       results(3,nchunk)= &_xyz3;
  363.     }
  364.   }
  365.   _xyz3= [];  /* free memory */
  366.  
  367.   /* collect the results of the chunking loop */
  368.   if (!ntotal) return [];
  369.   if (need_clist) clist= array(0, ntotal);
  370.   fs= array(0.0, 2,2,2,ntotal);
  371.   if (got_xyz) xyz= array(0.0, 3, dimsof(fs));
  372.   for (i=1,k=0 ; i<=nchunk ; ++i) {
  373.     j= k+1;
  374.     k+= numberof(*results(1,i));
  375.     if (need_clist) clist(j:k)= *results(1,i);
  376.     fs(,,,j:k)= *results(2,i);
  377.     if (!is_void(xyz)) xyz(,,,,j:k)= *results(3,i);
  378.   }
  379.   results= [];  /* free memory */
  380.   if (!got_xyz) xyz= xyz3(m3, clist);
  381.  
  382.   /* produce the lists of edge intersection points
  383.    * -- generate 2x2x3x(nsliced) array of edge mask values
  384.    * (mask non-zero if edge is cut by plane) */
  385.   below= fs<0.0;
  386.   mask= array(0n, 2,2,3,ntotal);
  387.   mask(-,,,1,)= abs(below(dif,,,));
  388.   mask(,-,,2,)= abs(below(,dif,,));
  389.   mask(,,-,3,)= abs(below(,,dif,));
  390.   list= where(mask);
  391.   edges= list-1;
  392.   cells= edges/12;      /* 0-origin momentarily */
  393.   edges= (edges%12)+1;
  394.   /* construct edge endpoint indices in fs, xyz arrays
  395.    * the numbers are the endpoint indices corresponding to
  396.    * the order of the 12 edges in the mask array */
  397.   lower= [1,3,5,7, 1,2,5,6, 1,2,3,4](edges) + 8*cells;
  398.   upper= [2,4,6,8, 3,4,7,8, 5,6,7,8](edges) + 8*cells;
  399.   /* restore the cells array to 1-origin */
  400.   cells+= 1;
  401.  
  402.   /* interpolate to find edge intersection points */
  403.   fsl= fs(lower)(-,);
  404.   fsu= fs(upper)(-,);
  405.   /* following denominator guaranteed non-zero */
  406.   xyz= (xyz(,lower)*fsu-xyz(,upper)*fsl)/(fsu-fsl);
  407.  
  408.   /* The xyz array is now the output xyzverts array,
  409.    * but for the order of the points within each cell.  */
  410.  
  411.   /* give each sliced cell a "pattern index" between 0 and 255
  412.    * (non-inclusive) representing the pattern of its 8 corners
  413.    * above and below the slicing plane */
  414.   pattern= below(*,)(+,) * (1<<indgen(0:7))(+);
  415.   if (slice3_stats) {
  416.     extern poly_patterns;
  417.     poly_patterns= histogram(pattern, top=254);
  418.   }
  419.   /* broadcast the cell's pattern onto each of its sliced edges */
  420.   pattern= pattern(-:1:12,)(list);
  421.  
  422.   /* To each pattern, there corresponds a permutation of the
  423.    * twelve edges so that they occur in the order in which the
  424.    * edges are to be connected.  Let each such permuation be
  425.    * stored as a list of integers from 0 to 11 such that sorting
  426.    * the integers into increasing order rearranges the edges at
  427.    * the corresponding indices into the correct order.  (The position
  428.    * of unsliced edges in the list is arbitrary as long as the sliced
  429.    * edges are in the proper order relative to each other.)
  430.    * Let these permutations be stored in a 12-by-254 array
  431.    * poly_permutations (see next comment for explanation of 48): */
  432.   pattern= poly_permutations(12*(pattern-1)+edges) + 48*cells;
  433.   order= sort(pattern);
  434.   xyz= xyz(,order);
  435.   edges= edges(order);
  436.   pattern= pattern(order);
  437.   /* cells(order) is same as cells by construction */
  438.  
  439.   /* There remains only the question of splitting the points in
  440.    * a single cell into multiple disjoint polygons.
  441.    * To do this, we need one more precomputed array: poly_splits
  442.    * should be another 12-by-254 array with values between 0 and 3
  443.    * 0 for each edge on the first part, 1 for each edge on the
  444.    * second part, and so on up to 3 for each edge on the fourth
  445.    * part.  The value on unsliced edges can be anything, say 0.
  446.    * With a little cleverness poly_splits can be combined with
  447.    * poly_permutations, by putting poly_permutations =
  448.    * poly_permutations(as described above) + 12*poly_splits
  449.    * (this doesn't change the ordering of poly_permutations).
  450.    * I assume this has been done here: */
  451.   pattern/= 12;
  452.   /* now pattern jumps by 4 between cells, smaller jumps within cells
  453.    * get the list of places where a new value begins, and form a
  454.    * new pattern with values that increment by 1 between each plateau */
  455.   pattern= pattern(dif);
  456.   list= grow([1], where(pattern)+1);
  457.   pattern= (pattern!=0)(cum) + 1;
  458.  
  459.   nverts= histogram(pattern);
  460.   xyzverts= xyz;
  461.  
  462.   /* finally, deal with any fcolor function */
  463.   if (is_void(fcolor)) return [];
  464.  
  465.   /* if some polys have been split, need to split clist as well */
  466.   if (numberof(list)>numberof(clist)) clist= clist(cells(list));
  467.   if (!nointerp) {
  468.     if (is_func(fcolor))
  469.       return fcolor(m3, clist, lower, upper, fsl(1,), fsu(1,), cells);
  470.     else
  471.       return getc3(fcolor, m3, clist, lower, upper, fsl(1,), fsu(1,), cells);
  472.   } else {
  473.     if (is_func(fcolor)) return fcolor(m3, clist);
  474.     else return getc3(fcolor, m3, clist);
  475.   }
  476. }
  477.  
  478. func _isosurface_slicer(m3, chunk)
  479. {
  480.   return getv3(iso_index, m3, chunk)-value;
  481. }
  482.  
  483. func _plane_slicer(m3, chunk)
  484. {
  485.   return xyz3(m3, chunk)(+,..)*normal(+) - projection;
  486. }
  487.  
  488. func to_corners3(list, ni, nj)
  489. /* DOCUMENT to_corners(list, ni, nj)
  490.      convert a LIST of cell indices in an (NI-1)-by-(NJ-1)-by-(nk-1)
  491.      logically rectangular grid of cells into the list of
  492.      2-by-2-by-2-by-numberof(LIST) cell corner indices in the
  493.      corresponding NI-by-NJ-by-nk list of vertices.
  494.  */
  495. {
  496.   ninj= ni*nj;
  497.   list-= 1;
  498.   ii= list/(ni-1);
  499.   list+= ii + ni*(ii/(nj-1));
  500.   return ([1,2]+[0,ni](-,)+[0,ninj](-,-,)) + list(-,-,-,);
  501. }
  502.  
  503. /* The iterator3 function combines three distinct operations:
  504.  * (1) If only the M3 argument is given, return the initial
  505.  *     chunk of the mesh.  The chunk will be no more than
  506.  *     chunk3_limit cells of the mesh.
  507.  * (2) If only M3 and CHUNK are given, return the next CHUNK,
  508.  *     or [] if there are no more chunks.
  509.  * (3) If M3, CHUNK, and CLIST are all specified, return the
  510.  *     absolute cell index list corresponding to the index list
  511.  *     CLIST of the cells in the CHUNK.
  512.  *     Do not increment the chunk in this case.
  513.  *
  514.  * The form of the CHUNK argument and return value for cases (1)
  515.  * and (2) is not specified, but it must be recognized by the
  516.  * xyz3 and getv3 functions which go along with this iterator3.
  517.  * (For case (3), CLIST and the return value are both ordinary
  518.  *  index lists.)
  519.  */
  520. func iterator3(m3, chunk, clist)
  521. {
  522.   return _car(_car(m3),4)(m3, chunk, clist);
  523. }
  524.  
  525. /* biggest temporary is 3 doubles times this,
  526.  * perhaps 4 or 5 doubles times this is most at one time */
  527. chunk3_limit= 10000;
  528.  
  529. func iterator3_rect(m3, chunk, clist)
  530. {
  531.   if (is_void(chunk)) {
  532.     dims= _car(_car(m3,2));  /* [ni,nj,nk] cell dimensions */
  533.     ni= dims(1);
  534.     nj= dims(2);
  535.     nk= dims(3);
  536.     ninj= ni*nj;
  537.     if (chunk3_limit <= ni) {
  538.       /* stuck with 1D chunks */
  539.       ci= (ni-1)/chunk3_limit + 1;
  540.       cj= ck= 0;
  541.     } else if (chunk3_limit <= ninj) {
  542.       /* 2D chunks */
  543.       ci= ck= 0;
  544.       cj= (ninj-1)/chunk3_limit + 1;
  545.     } else {
  546.       /* 3D chunks */
  547.       ci= cj= 0;
  548.       ck= (ninj*nk-1)/chunk3_limit + 1;
  549.     }
  550.     chunk= [[ci==0,cj==0,ck==0],
  551.         [ni*((cj+ck)!=0),nj*(ck!=0)+(ci!=0),!ck],
  552.         [ci,cj,ck],[ni,nj,nk]];
  553.  
  554.   } else {
  555.     chunk= *chunk;
  556.     ni= chunk(1,4);  nj= chunk(2,4);  nk= chunk(3,4);
  557.     ninj= ni*nj;
  558.  
  559.     if (!is_void(clist)) {
  560.       /* add offset for this chunk to clist and return */
  561.       return [1,ni,ninj](+)*(chunk(,1)-1)(+) + clist;
  562.     }
  563.   }
  564.  
  565.   /* increment to next chunk */
  566.   xi= chunk(1,2);  xj= chunk(2,2);  xk= chunk(3,2);
  567.  
  568.   if ((np= chunk(1,3))) {
  569.     /* 1D chunks */
  570.     if (xi==ni) {
  571.       if (xj==nj) {
  572.     if (xk==nk) return [];
  573.     xk+= 1;
  574.     xj= 1;
  575.       } else {
  576.     xj+= 1;
  577.       }
  578.       xi= 0;
  579.     }
  580.     ci= xi+1;
  581.     step= ni/np;
  582.     frst= ni%np;   /* first frst steps are step+1 */
  583.     if (xi < (step+1)*frst) step+= 1;
  584.     xi+= step;
  585.     chunk(,1)= [ci,xj,xk];
  586.     chunk(,2)= [xi,xj,xk];
  587.  
  588.   } else if ((np= chunk(2,3))) {
  589.     if (xj==nj) {
  590.       if (xk==nk) return [];
  591.       xk+= 1;
  592.       xj= 0;
  593.     }
  594.     cj= xj+1;
  595.     step= nj/np;
  596.     frst= nj%np;   /* first frst steps are step+1 */
  597.     if (xj < (step+1)*frst) step+= 1;
  598.     xj+= step;
  599.     chunk(2:3,1)= [cj,xk];
  600.     chunk(2:3,2)= [xj,xk];
  601.  
  602.   } else {
  603.     if (xk==nk) return [];
  604.     ck= xk+1;
  605.     np= chunk(3,3);
  606.     step= nk/np;
  607.     frst= nk%np;   /* first frst steps are step+1 */
  608.     if (xk < (step+1)*frst) step+= 1;
  609.     xk+= step;
  610.     chunk(3,1)= ck;
  611.     chunk(3,2)= xk;
  612.   }
  613.  
  614.   /* return pointer to make chunk easy to recognize for xyz3, getv3 */
  615.   return &chunk;
  616. }
  617.  
  618. func xyz3(m3, chunk)
  619. /* DOCUMENT xyz3(m3, chunk)
  620.  
  621.      return vertex coordinates for CHUNK of 3D mesh M3.  The CHUNK
  622.      may be a list of cell indices, in which case xyz3 returns a
  623.      3x2x2x2x(dimsof(CHUNK)) list of vertex coordinates.  CHUNK may
  624.      also be a mesh-specific data structure used in the slice3
  625.      routine, in which case xyz3 may return a 3x(ni)x(nj)x(nk)
  626.      array of vertex coordinates.  For meshes which are logically
  627.      rectangular or consist of several rectangular patches, this
  628.      is up to 8 times less data, with a concomitant performance
  629.      advantage.  Use xyz3 when writing slicing functions or coloring
  630.      functions for slice3.
  631.  
  632.    SEE ALSO: slice3, mesh3, getv3, getc3
  633.  */
  634. {
  635.   xyz= _car(_car(m3),1)(m3, chunk);
  636.   extern _xyz3;
  637.   /* the cost of the following copy operation could be saved
  638.    * by making _xyz3 a pointer, but then we would be relying on
  639.    * the fslice caller of xyz3 not using the returned array for
  640.    * scratch space - leave it as is for now */
  641.   if (_xyz3_save) _xyz3= xyz;
  642.   return xyz;
  643. }
  644.  
  645. func xyz3_rect(m3, chunk)
  646. {
  647.   m3= _car(m3,2);
  648.   if (structof(chunk)==pointer) {
  649.     c= *chunk;
  650.     return _car(m3,2)(,c(1,1):1+c(1,2),c(2,1):1+c(2,2),c(3,1):1+c(3,2));
  651.   } else {
  652.     dims= _car(m3);
  653.     return _car(m3,2)(,to_corners3(chunk,dims(1)+1,dims(2)+1));
  654.   }
  655. }
  656.  
  657. func xyz3_unif(m3, chunk)
  658. {
  659.   m3= _car(m3,2);
  660.   n= _car(m3,2);
  661.   dxdydz= n(,1);
  662.   x0y0z0= n(,2);
  663.   if (structof(chunk)==pointer) {
  664.     c= *chunk;
  665.     i= c(,1)-1;
  666.     n= c(,2)+1-i;
  667.     xyz= array(x0y0z0+i*dxdydz, n(1), n(2), n(3));
  668.     xyz(1,..)+= span(0.,dxdydz(1),n(1));
  669.     xyz(2,..)+= span(0.,dxdydz(2),n(2))(-,);
  670.     xyz(3,..)+= span(0.,dxdydz(3),n(3))(-,-,);
  671.     return xyz;
  672.   } else {
  673.     dims= _car(m3);
  674.     ni= dims(1);
  675.     nj= dims(2);
  676.     ninj= ni*nj;
  677.     ijk= (chunk-1)(-:1:3,-,-,-,..);
  678.     ijk(3,..)/= ninj;
  679.     ijk(2,..)%= ninj;
  680.     ijk(2,..)/= ni;
  681.     ijk(1,..)%= ni;
  682.     ijk+= [[[[0,0,0],[1,0,0]],[[0,1,0],[1,1,0]]],
  683.        [[[0,0,1],[1,0,1]],[[0,1,1],[1,1,1]]]];
  684.     return x0y0z0 + ijk*dxdydz;
  685.   }
  686. }
  687.  
  688. func getv3(i, m3, chunk)
  689. /* DOCUMENT getv3(i, m3, chunk)
  690.  
  691.      return vertex values of the Ith function attached to 3D mesh M3
  692.      for cells in the specified CHUNK.  The CHUNK may be a list of
  693.      cell indices, in which case getv3 returns a 2x2x2x(dimsof(CHUNK))
  694.      list of vertex coordinates.  CHUNK may also be a mesh-specific data
  695.      structure used in the slice3 routine, in which case getv3 may
  696.      return a (ni)x(nj)x(nk) array of vertex values.  For meshes which
  697.      are logically rectangular or consist of several rectangular
  698.      patches, this is up to 8 times less data, with a concomitant
  699.      performance advantage.  Use getv3 when writing slicing functions
  700.      for slice3.
  701.  
  702.    SEE ALSO: slice3, mesh3, getc3, xyz3
  703.  */
  704. {
  705.   return _car(_car(m3),2)(i, m3, chunk);
  706. }
  707.  
  708. func getv3_rect(i, m3, chunk)
  709. {
  710.   fi= _cdr(m3,2);
  711.   m3= _car(m3,2);
  712.   if (i<1 || i>_len(fi)) error, "no such mesh function as F"+pr1(i);
  713.   dims= _car(m3);
  714.   if (allof(dimsof(_car(fi,i))(2:4)==dims)) {
  715.     error, "mesh function F"+pr1(i)+" is not vertex-centered";
  716.   }
  717.   if (structof(chunk)==pointer) {
  718.     c= *chunk;
  719.     return _car(fi,i)(c(1,1):1+c(1,2),c(2,1):1+c(2,2),c(3,1):1+c(3,2));
  720.   } else {
  721.     return _car(fi,i)(to_corners3(chunk,dims(1)+1,dims(2)+1));
  722.   }
  723. }
  724.  
  725. func getc3(i, m3, chunk, l, u, fsl, fsu, cells)
  726. /* DOCUMENT getc3(i, m3, chunk)
  727.          or getc3(i, m3, clist, l, u, fsl, fsu, cells)
  728.  
  729.      return cell values of the Ith function attached to 3D mesh M3
  730.      for cells in the specified CHUNK.  The CHUNK may be a list of
  731.      cell indices, in which case getc3 returns a (dimsof(CHUNK))
  732.      list of vertex coordinates.  CHUNK may also be a mesh-specific data
  733.      structure used in the slice3 routine, in which case getc3 may
  734.      return a (ni)x(nj)x(nk) array of vertex values.  There is no
  735.      savings in the amount of data for such a CHUNK, but the gather
  736.      operation is cheaper than a general list of cell indices.
  737.      Use getc3 when writing colorng functions for slice3.
  738.  
  739.      If CHUNK is a CLIST, the additional arguments L, U, FSL, and FSU
  740.      are vertex index lists which override the CLIST if the Ith attached
  741.      function is defined on mesh vertices.  L and U are index lists into
  742.      the 2x2x2x(dimsof(CLIST)) vertex value array, say vva, and FSL
  743.      and FSU are corresponding interpolation coefficients; the zone
  744.      centered value is computed as a weighted average of involving these
  745.      coefficients.  The CELLS argument is required by histogram to do
  746.      the averaging.  See the source code for details.
  747.      By default, this conversion (if necessary) is done by averaging
  748.      the eight vertex-centered values.
  749.  
  750.    SEE ALSO: slice3, mesh3, getv3, xyz3
  751.  */
  752. {
  753.   return _car(_car(m3),3)(i, m3, chunk, l, u, fsl, fsu, cells);
  754. }
  755.  
  756. func getc3_rect(i, m3, chunk, l, u, fsl, fsu, cells)
  757. {
  758.   fi= _cdr(m3,2);
  759.   m3= _car(m3,2);
  760.   if (i<1 || i>_len(fi)) error, "no such mesh function as F"+pr1(i);
  761.   dims= _car(m3);
  762.   if (allof(dimsof(_car(fi,i))(2:4)==dims)) {
  763.     if (structof(chunk)==pointer) {
  764.       c= *chunk;
  765.       return _car(fi,i)(c(1,1):c(1,2),c(2,1):c(2,2),c(3,1):c(3,2));
  766.     } else {
  767.       return _car(fi,i)(chunk);
  768.     }
  769.   } else {
  770.     if (structof(chunk)==pointer) {
  771.       c= *chunk;
  772.       return _car(fi,i)(zcen:c(1,1):1+c(1,2),zcen:c(2,1):1+c(2,2),
  773.                 zcen:c(3,1):1+c(3,2));
  774.     } else {
  775.       corners= _car(fi,i)(to_corners3(chunk,dims(1)+1,dims(2)+1));
  776.       if (is_void(l)) {
  777.     return 0.125*corners(sum:1:8,1,1,..);
  778.       } else {
  779.     /* interpolate corner values to get edge values */
  780.     corners= (corners(l)*fsu-corners(u)*fsl)/(fsu-fsl);
  781.     /* average edge values (vertex values of polys) on each poly */
  782.     return histogram(cells,corners)/histogram(cells);
  783.       }
  784.     }
  785.   }
  786. }
  787.  
  788. /* ------------------------------------------------------------------------ */
  789.  
  790. /* There are 254 possible combinations of 8 vertices above and below
  791.  * the slicing plane (not counting all above or all below).  */
  792. func _construct3(void)
  793. {
  794.   /* construct the edge list for each of the 254 possibilities */
  795.   i= indgen(254);
  796.   below= transpose([[[i&1,i&2],[i&4,i&8]],[[i&16,i&32],[i&64,i&128]]],0);
  797.   below= (below!=0);
  798.   mask= array(0n, 2,2,3,254);
  799.   mask(-,,,1,)= abs(below(dif,,,));
  800.   mask(,-,,2,)= abs(below(,dif,,));
  801.   mask(,,-,3,)= abs(below(,,dif,));
  802.  
  803.   /* walking the edges requires some connectivity information */
  804.   start_face= [3,4,3,4, 1,2,1,2, 1,2,1,2];
  805.   face_edges= [[5,11,7,9], [6,10,8,12],
  806.            [1,9,3,10], [2,12,4,11],
  807.            [1,6,2,5],  [3,7,4,8]];
  808.   edge_faces= [[3,5],[4,5],[3,6],[4,6],[1,5],[2,5],
  809.            [1,6],[2,6],[1,3],[2,3],[1,4],[2,4]];
  810.  
  811.   permute= array('\0', 12, 254);
  812.   for (i=1 ; i<=254 ; ++i) permute(,i)= _walk3(mask(*,i));
  813.   return permute;
  814. }
  815.  
  816. func _walk3(mask)
  817. {
  818.   permute= splits= array(0, 12);
  819.   split= 0;
  820.   list= where(mask);
  821.   edge= min(list);
  822.   face= start_face(edge);
  823.  
  824.   for (i=0 ; i<numberof(list)-1 ; ++i) {
  825.     /* record this edge */
  826.     permute(edge)= i;
  827.     splits(edge)= split;
  828.     mask(edge)= 0;
  829.  
  830.     /* advance to next edge */
  831.     try= face_edges(,face);
  832.     now= abs(try-edge)(mnx);
  833.     /* test opposite edge first (make X, never // or \\) */
  834.     edge= try((now+1)%4+1);
  835.     if (!mask(edge)) {
  836.       /* otherwise one or the other opposite edges must be set */
  837.       edge= try(now%4+1);
  838.       if (!mask(edge)) {
  839.     edge= try((now+2)%4+1);
  840.     if (!mask(edge)) {
  841.       /* we must have closed a loop */
  842.       split+= 1;
  843.       edge= min(where(mask));
  844.     }
  845.       }
  846.     }
  847.     try= edge_faces(,edge);
  848.     now= try(1);
  849.     face= face==now? try(2) : now;
  850.   }
  851.   /* last pass through loop */
  852.   permute(edge)= i;
  853.   splits(edge)= split;
  854.   mask(edge)= 0;
  855.  
  856.   return permute+12*splits;
  857. }
  858.  
  859. poly_permutations= _construct3();
  860.  
  861. /* ------------------------------------------------------------------------ */
  862.  
  863. func slice2x(plane, &nverts, &xyzverts, &values, &nvertb, &xyzvertb, &valueb)
  864. /* DOCUMENT slice2, plane, nverts, values, xyzverts
  865.  
  866.      Slice a polygon list, retaining only those polygons or
  867.      parts of polygons on the positive side of PLANE, that is,
  868.      the side where xyz(+)*PLANE(+:1:3)-PLANE(4) > 0.0.
  869.      The NVERTS, VALUES, and XYZVERTS arrays serve as both
  870.      input and output, and have the meanings of the return
  871.      values from the slice3 function.
  872.  
  873.    SEE ALSO: slice2, slice2_precision
  874.  */
  875. {
  876.   _slice2x= 1;
  877.   return slice2(plane, nverts, xyzverts, values, nvertb, xyzvertb, valueb);
  878. }
  879.  
  880. func slice2(plane, &nverts, &xyzverts, &values, &nvertb, &xyzvertb, &valueb)
  881. /* DOCUMENT slice2, plane, nverts, xyzverts
  882.          or slice2, plane, nverts, xyzverts, values
  883.  
  884.      Slice a polygon list, retaining only those polygons or
  885.      parts of polygons on the positive side of PLANE, that is,
  886.      the side where xyz(+)*PLANE(+:1:3)-PLANE(4) > 0.0.
  887.      The NVERTS, XYZVERTS, and VALUES arrays serve as both
  888.      input and output, and have the meanings of the return
  889.      values from the slice3 function.  It is legal to omit the
  890.      VALUES argument (e.g.- if there is no fcolor function).
  891.  
  892.      In order to plot two intersecting slices, one could
  893.      slice (for example) the horizontal plane twice (slice2x) -
  894.      first with the plane of the vertical slice, then with minus
  895.      that same plane.  Then, plot first the back part of the
  896.      slice, then the vertical slice, then the front part of the
  897.      horizontal slice.  Of course, the vertical plane could
  898.      be the one to be sliced, and "back" and "front" vary
  899.      depending on the view point, but the general idea always
  900.      works.
  901.  
  902.    SEE ALSO: slice3, plane3, slice2x, slice2_precision
  903.  */
  904. {
  905.   have_values= !is_void(values);
  906.  
  907.   /* get the list of indices into nverts (or values) for each of
  908.    * the points in xyzverts */
  909.   ndxs= histogram(1+nverts(psum))(1:-1);
  910.   ndxs(1)+= 1;
  911.   ndxs= ndxs(psum);
  912.  
  913.   /* form dot products of all the points with the given plane */
  914.   dp= xyzverts(+,)*plane(+:1:3) - plane(4);
  915.  
  916.   /* separate into lists of unclipped and partially clipped polys */
  917.   if (!slice2_precision) {
  918.     /* if precision is not set, slice exactly at dp==0.0, with
  919.      * any points exactly at dp==0.0 treated as if they had dp>0.0 */
  920.     keep= (dp>=0.0);
  921.   } else {
  922.     /* if precision is set, polygons are clipped to +-precision,
  923.      * so that any poly crossing +precision is clipped to dp>=+precision,
  924.      * any poly crossing -precision is clipped to dp<=-precision, and
  925.      * any poly lying entirely between +-precision is discarded entirely */
  926.     keep= (dp>=slice2_precision);
  927.   }
  928.   nkeep= long(histogram(ndxs, keep));
  929.   mask0= (nkeep==nverts);
  930.   mask1= (nkeep!=0 & !mask0);
  931.   list1= where(mask1);
  932.   if (numberof(list1)) {
  933.     nvertc= nverts(list1);
  934.     if (have_values) valuec= values(list1);
  935.     list= where(mask1(ndxs));
  936.     xyzc= xyzverts(, list);
  937.   }
  938.   if (_slice2x) {
  939.     if (!slice2_precision) {
  940.       mask2= !nkeep;
  941.       nvertc0= nvertc;
  942.       valuec0= valuec;
  943.       xyzc0= xyzc;
  944.     } else {
  945.       keep2= (dp>-slice2_precision);
  946.       nkeep2= long(histogram(ndxs, keep2));
  947.       mask2= (nkeep!=0 & nkeep<nverts);
  948.       list2= where(mask2);
  949.       if (numberof(list2)) {
  950.     nvertc0= nverts(list2);
  951.     if (have_values) valuec0= values(list2);
  952.     listc= where(mask2(ndxs));
  953.     xyzc0= xyzverts(, listc);
  954.       }
  955.       mask2= !nkeep2;
  956.     }
  957.     list2= where(mask2);
  958.     if (numberof(list2)) {
  959.       nvertb= nverts(list2);
  960.       if (have_values) valueb= values(list2);
  961.       xyzvertb= xyzverts(, where(mask2(ndxs)));
  962.     } else {
  963.       nvertb= valueb= xyzvertb= [];
  964.     }
  965.   }
  966.   list0= where(mask0);
  967.   if (numberof(list0)<numberof(nverts)) {
  968.     nverts= nverts(list0);
  969.     if (have_values) values= values(list0);
  970.     xyzverts= xyzverts(, where(mask0(ndxs)));
  971.   }
  972.  
  973.   /* done if no partially clipped polys */
  974.   if (!numberof(list1) && !numberof(listc)) return;
  975.   if (!numberof(list1)) goto skip;
  976.  
  977.   /* get dot products and keep list for the clipped polys */
  978.   dp= dp(list);
  979.   if (slice2_precision) dp-= slice2_precision;
  980.   keep= (dp>=0.0);
  981.  
  982.   /* get the indices of the first and last points in each clipped poly */
  983.   last= nvertc(psum);
  984.   frst= last - nvertc + 1;
  985.  
  986.   /* get indices of previous and next points in cyclic order */
  987.   prev= next= indgen(0:numberof(keep)-1);
  988.   prev(frst)= last;
  989.   next(prev)= indgen(numberof(keep));
  990.  
  991.   _slice2_part;
  992.  
  993.   grow, nverts, nvertc;
  994.   if (have_values) grow, values, valuec;
  995.   grow, xyzverts, xyzc;
  996.  
  997.   if (_slice2x) {
  998.   skip:
  999.     nvertc= nvertc0;
  1000.     valuec= valuec0;
  1001.     xyzc= xyzc0;
  1002.     if (!slice2_precision) {
  1003.       keep= !keep;
  1004.     } else {
  1005.       dp= dp(listc)+slice2_precision;
  1006.       keep= (dp>=0.0);
  1007.     }
  1008.  
  1009.     _slice2_part;
  1010.  
  1011.     grow, nvertb, nvertc;
  1012.     if (have_values) grow, valueb, valuec;
  1013.     grow, xyzvertb, xyzc;
  1014.   }
  1015. }
  1016.  
  1017. local slice2_precision;
  1018. /* DOCUMENT slice2_precision= precision
  1019.      Controls how slice2 (or slice2x) handles points very close to
  1020.      the slicing plane.  PRECISION should be a positive number or zero.
  1021.      Zero PRECISION means to clip exactly to the plane, with points
  1022.      exactly on the plane acting as if they were slightly on the side
  1023.      the normal points toward.  Positive PRECISION means that edges
  1024.      are clipped to parallel planes a distance PRECISION on either
  1025.      side of the given plane.  (Polygons lying entirely between these
  1026.      planes are completely discarded.)
  1027.  
  1028.      Default value is 0.0.
  1029.  
  1030.    SEE ALSO: slice2, slice2x
  1031.  */
  1032. slice2_precision= 0.0;
  1033.  
  1034. func _slice2_part
  1035. {
  1036.   extern nvertc, xyzc;
  1037.  
  1038.   /* find the points where any edges of the polys cut the plane */
  1039.   mask0= (!keep) & keep(next);   /* off-->on */
  1040.   list0= where(mask0);
  1041.   if (numberof(list0)) {
  1042.     list= next(list0);
  1043.     dpl= dp(list0)(-,);
  1044.     dpu= dp(list)(-,);
  1045.     xyz0= (xyzc(,list0)*dpu-xyzc(,list)*dpl)/(dpu-dpl);
  1046.   }
  1047.   mask1= (!keep) & keep(prev);   /* on-->off */
  1048.   list1= where(mask1);
  1049.   if (numberof(list1)) {
  1050.     list= prev(list1);
  1051.     dpl= dp(list1)(-,);
  1052.     dpu= dp(list)(-,);
  1053.     xyz1= (xyzc(,list1)*dpu-xyzc(,list)*dpl)/(dpu-dpl);
  1054.   }
  1055.  
  1056.   /* form an index list xold which gives the indices in the original
  1057.    * xyzc list at the places in the new, sliced xyzc list */
  1058.   mask= keep+mask0+mask1;  /* 0, 1, or 2 */
  1059.   list= mask(psum);  /* index values in new list */
  1060.   xold= array(0, list(0));
  1061.   mlist= where(mask);
  1062.   xold(list(mlist))= mlist;
  1063.   dups= where(mask==2);
  1064.   if (numberof(dups)) xold(list(dups)-1)= dups;
  1065.  
  1066.   /* form the new, sliced xyzc vertex list */
  1067.   xyzc= xyzc(,xold);
  1068.   if (numberof(list0)) xyzc(,list(list0))= xyz0;
  1069.   if (numberof(list1)) xyzc(,list(list1)-mask(list1)+1)= xyz1;
  1070.  
  1071.   /* get the list of indices into nvertc (or valuec) for each of
  1072.    * the points in xyzc */
  1073.   ndxs= histogram(1+last)(1:-1);
  1074.   ndxs(1)+= 1;
  1075.   ndxs= ndxs(psum);
  1076.   /* compute the new number of vertices */
  1077.   nvertc= histogram(ndxs(xold));
  1078. }
  1079.  
  1080. /* ------------------------------------------------------------------------ */
  1081.  
  1082. func pl3surf(nverts, xyzverts, values, cmin=, cmax=)
  1083. /* DOCUMENT pl3surf, nverts, xyzverts
  1084.          or pl3surf, nverts, xyzverts, values
  1085.  
  1086.      Perform simple 3D rendering of an object created by slice3
  1087.      (possibly followed by slice2).  NVERTS and XYZVERTS are polygon
  1088.      lists as returned by slice3, so XYZVERTS is 3-by-sum(NVERTS),
  1089.      where NVERTS is a list of the number of vertices in each polygon.
  1090.      If present, the VALUES should have the same length as NVERTS;
  1091.      they are used to color the polygon.  If VALUES is not specified,
  1092.      the 3D lighting calculation set up using the light3 function
  1093.      will be carried out.  Keywords cmin= and cmax= as for plf, pli,
  1094.      or plfp are also accepted.  (If you do not supply VALUES, you
  1095.      probably want to use the ambient= keyword to light3 instead of
  1096.      cmin= here, but cmax= may still be useful.)
  1097.  
  1098.    SEE ALSO: pl3tree, slice3, slice2, rot3, light3
  1099.  */
  1100. {
  1101.   require, "pl3d.i";
  1102.   if (_draw3) {
  1103.     list= nverts;
  1104.     nverts= _nxt(list);
  1105.     xyzverts= _nxt(list);
  1106.     values= _nxt(list);
  1107.     cmin= _nxt(list);
  1108.     cmax= _nxt(list);
  1109.  
  1110.     local x, y, z, list, vlist;
  1111.     get3_xy, xyzverts, x, y, z, 1;
  1112.     if (is_void(values)) {
  1113.       xyzverts(1,..)= x;
  1114.       xyzverts(2,..)= y;
  1115.       xyzverts(3,..)= z;
  1116.       values= get3_light(xyzverts, nverts);
  1117.     }
  1118.     sort3d, z, nverts, list, vlist;
  1119.     nverts= nverts(list);
  1120.     values= values(list);
  1121.     x= x(vlist);
  1122.     y= y(vlist);
  1123.  
  1124.     plfp, values,y,x,nverts, cmin=cmin,cmax=cmax, legend=string(0);
  1125.     return;
  1126.   }
  1127.  
  1128.   nverts+= 0;
  1129.   xyzverts+= 0.0;
  1130.  
  1131.   if (numberof(xyzverts)!=3*sum(nverts) || anyof(nverts<3) ||
  1132.       structof(nverts)!=long)
  1133.     error, "illegal or inconsistent polygon list";
  1134.   if (!is_void(values) && numberof(values)!=numberof(nverts))
  1135.     error, "illegal or inconsistent polygon color values";
  1136.  
  1137.   if (!is_void(values)) values+= 0.0;
  1138.  
  1139.   clear3;
  1140.   set3_object, pl3surf, _lst(nverts,xyzverts,values,cmin,cmax);
  1141. }
  1142.  
  1143. /* ------------------------------------------------------------------------ */
  1144.  
  1145. func pl3tree(nverts, xyzverts, values, plane, cmin=, cmax=)
  1146. /* DOCUMENT pl3tree, nverts, xyzverts
  1147.          or pl3tree, nverts, xyzverts, values, plane
  1148.  
  1149.      Add the polygon list specified by NVERTS (number of vertices in
  1150.      each polygon) and XYZVERTS (3-by-sum(NVERTS) vertex coordinates)
  1151.      to the currently displayed b-tree.  If VALUES is specified, it
  1152.      must have the same dimension as NVERTS, and represents the color
  1153.      of each polygon.  If VALUES is not specified, the polygons
  1154.      are assumed to form an isosurface which will be shaded by the
  1155.      current 3D lighting model; the isosurfaces are at the leaves of
  1156.      the b-tree, sliced by all of the planes.  If PLANE is specified,
  1157.      the XYZVERTS must all lie in that plane, and that plane becomes
  1158.      a new slicing plane in the b-tree.  
  1159.  
  1160.      Each leaf of the b-tree consists of a set of sliced isosurfaces.
  1161.      A node of the b-tree consists of some polygons in one of the
  1162.      planes, a b-tree or leaf entirely on one side of that plane, and
  1163.      a b-tree or leaf on the other side.  The first plane you add
  1164.      becomes the root node, slicing any existing leaf in half.  When
  1165.      you add an isosurface, it propagates down the tree, getting
  1166.      sliced at each node, until its pieces reach the existing leaves,
  1167.      to which they are added.  When you add a plane, it also propagates
  1168.      down the tree, getting sliced at each node, until its pieces
  1169.      reach the leaves, which it slices, becoming the nodes closest to
  1170.      the leaves.
  1171.  
  1172.      This structure is relatively easy to plot, since from any
  1173.      viewpoint, a node can always be plotted in the order from one
  1174.      side, then the plane, then the other side.
  1175.  
  1176.      This routine assumes a "split palette"; the colors for the
  1177.      VALUES will be scaled to fit from color 0 to color 99, while
  1178.      the colors from the shading calculation will be scaled to fit
  1179.      from color 100 to color 199.  (If VALUES is specified as a char
  1180.      array, however, it will be used without scaling.)
  1181.      You may specifiy a cmin= or cmax= keyword to affect the
  1182.      scaling; cmin is ignored if VALUES is not specified (use the
  1183.      ambient= keyword from light3 for that case).
  1184.  
  1185.    SEE ALSO: pl3surf, slice3, slice2, rot3, light3, split_palette
  1186.  */
  1187. {
  1188.   require, "pl3d.i";
  1189.   if (_draw3) {
  1190.     /* avoid overhead of local variables for _pl3tree and _pl3leaf */
  1191.     local x,y,z;
  1192.     local nverts, xyzverts, values, list, vlist, cmin, cmax;
  1193.     local nv, vv, xx, yy, zz, cmin, cmax;
  1194.     _pl3tree, nverts;
  1195.     return;
  1196.   }
  1197.  
  1198.   nverts= 0+nverts(*);
  1199.   xyzverts= double(xyzverts(,*));
  1200.   if (!is_void(values)) values= double(values(*));
  1201.   if (!is_void(plane)) plane= double(plane);
  1202.  
  1203.   if (numberof(xyzverts)!=3*sum(nverts) || anyof(nverts<3) ||
  1204.       structof(nverts)!=long)
  1205.     error, "illegal or inconsistent polygon list";
  1206.   if (!is_void(values) && numberof(values)!=numberof(nverts))
  1207.     error, "illegal or inconsistent polygon color values";
  1208.   if (!is_void(plane) && (dimsof(plane)(1)!=1 || numberof(plane)!=4))
  1209.     error, "illegal plane format, try plane3 function";
  1210.  
  1211.   leaf= _lst(_lst(nverts, xyzverts, values, cmin, cmax));
  1212.  
  1213.   /* retrieve current b-tree (if any) from 3D display list */
  1214.   tree= _cdr(_draw3_list, _draw3_n);
  1215.   if (is_void(tree) || _car(tree)!=pl3tree) {
  1216.     tree= _lst(plane, [], leaf, []);
  1217.   } else {
  1218.     tree= _car(tree,2);
  1219.     _pl3tree_add, leaf, plane, tree;
  1220.   }
  1221.  
  1222.   clear3;
  1223.   set3_object, pl3tree, tree;
  1224. }
  1225.  
  1226. func split_palette(name)
  1227. /* DOCUMENT split_palette
  1228.          or split_palette, "palette_name.gp"
  1229.      split the current palette or the specified palette into two
  1230.      parts; colors 0 to 99 will be a compressed version of the
  1231.      original, while colors 100 to 199 will be a gray scale.
  1232.    SEE ALSO: pl3tree, split_bytscl
  1233.  */
  1234. {
  1235.   if (!is_void(name)) palette, name;
  1236.   local r,g,b;
  1237.   palette,query=1, r,g,b;
  1238.   n= numberof(r);
  1239.   if (n<100) {
  1240.     palette, "earth.gp";
  1241.     palette,query=1, r,g,b;
  1242.     n= numberof(r);
  1243.   }
  1244.   r= char(interp(r,indgen(n),span(1,n,100)));
  1245.   g= char(interp(g,indgen(n),span(1,n,100)));
  1246.   b= char(interp(b,indgen(n),span(1,n,100)));
  1247.   grow, r, char(span(0,255,100));
  1248.   grow, g, char(span(0,255,100));
  1249.   grow, b, char(span(0,255,100));
  1250.   palette, r,g,b;
  1251. }
  1252.  
  1253. func split_bytscl(x, upper, cmin=,cmax=)
  1254. /* DOCUMENT split_bytscl(x, 0)
  1255.          or split_bytscl(x, 1)
  1256.      as bytscl function, but scale to the lower half of a split
  1257.      palette (0-99, normally the color scale) if the second parameter
  1258.      is zero or nil, or the upper half (100-199, normally the gray
  1259.      scale) if the second parameter is non-zero.
  1260.    SEE ALSO: split_palette
  1261.  */
  1262. {
  1263.   x= bytscl(x,cmin=cmin,cmax=cmax,top=99);
  1264.   if (upper) x+= char(100);
  1265.   return x;
  1266. }
  1267.  
  1268. func _pl3tree(tree)
  1269. {
  1270.   /* avoid overhead of local variables for _pl3tree and _pl3leaf */
  1271.   extern x,y,z;
  1272.  
  1273.   /* tree is a 4-element list like this:
  1274.    * _lst(plane, back_tree, inplane_leaf, front_tree)
  1275.    * plane= _car(tree)  is void if this is just a leaf
  1276.    *                    in which case, only inplane_leaf is not void
  1277.    * back_tree= _car(tree,2)    is the part behind plane
  1278.    * inplane_leaf= _car(tree,3) is the part in the plane itself
  1279.    * front_tree= _car(tree,4)   is the part in front of plane
  1280.    */
  1281.   if (is_void(tree)) return;
  1282.   if (is_void(_car(tree))) {
  1283.     /* only the leaf is non-nil (but not a plane) */
  1284.     _pl3leaf, _car(tree,3), 1;
  1285.     return;
  1286.   }
  1287.  
  1288.   /* apply the 3D coordinate transform to two points along the
  1289.    * normal of the splitting plane to judge which is in front */
  1290.   get3_xy, [[0.,0.,0.],_car(tree)(1:3)], x,y,z, 1;
  1291.  
  1292.   /* plot the parts in order toward the current viewpoint */
  1293.   if (z(2) >= z(1)) {
  1294.     _pl3tree, _car(tree,2);
  1295.     _pl3leaf, _car(tree,3), 0;
  1296.     _pl3tree, _car(tree,4);
  1297.   } else {
  1298.     _pl3tree, _car(tree,4);
  1299.     _pl3leaf, _car(tree,3), 0;
  1300.     _pl3tree, _car(tree,2);
  1301.   }
  1302. }
  1303.  
  1304. func _pl3leaf(leaf, not_plane)
  1305. {
  1306.   /* avoid overhead of local variables for _pl3tree and _pl3leaf */
  1307.   extern x,y,z;
  1308.   extern nverts, xyzverts, values, list, vlist;
  1309.  
  1310.   /* count number of polys, number of vertices */
  1311.   nverts= xyzverts= 0;
  1312.   _map, _pl3tree_count, leaf;
  1313.  
  1314.   /* accumulate polys and vertices into a single polygon list */
  1315.   values= array(char, nverts);
  1316.   nverts= array(0, nverts);
  1317.   x= y= array(0.0, xyzverts);
  1318.   if (not_plane) z= x;
  1319.   else z= [];
  1320.   list= vlist= 1;
  1321.   _map, _pl3tree_accum, leaf;
  1322.  
  1323.   /* sort the single polygon list */
  1324.   if (not_plane) {
  1325.     sort3d, z, nverts, list, vlist;
  1326.     nverts= nverts(list);
  1327.     values= values(list);
  1328.     x= x(vlist);
  1329.     y= y(vlist);
  1330.   }
  1331.  
  1332.   plfp, values,y,x,nverts, legend=string(0);
  1333. }
  1334.  
  1335. func _pl3tree_count(item)
  1336. {
  1337.   nverts+= numberof(_nxt(item));
  1338.   xyzverts+= numberof(_nxt(item))/3;
  1339. }
  1340.  
  1341. func _pl3tree_accum(item)
  1342. {
  1343.   /* avoid overhead of local variables for _pl3tree and _pl3leaf */
  1344.   extern x,y,z;
  1345.   extern nverts, xyzverts, values, list, vlist;
  1346.   extern nv, vv, xx, yy, zz, cmin, cmax;
  1347.  
  1348.   nv= _nxt(item);
  1349.   xyzverts= _nxt(item);
  1350.   vv= _nxt(item);
  1351.   cmin= _nxt(item);
  1352.   cmax= _nxt(item);
  1353.  
  1354.   if (is_void(vv)) {
  1355.     /* this is an isosurface to be shaded */
  1356.     get3_xy, xyzverts, xx, yy, zz, 1;
  1357.     xyzverts(1,..)= xx;
  1358.     xyzverts(2,..)= yy;
  1359.     xyzverts(3,..)= zz;
  1360.     vv= get3_light(xyzverts, nv);
  1361.     vv= split_bytscl(vv, 1, cmin=0.0, cmax=cmax);
  1362.   } else {
  1363.     /* this is to be pseudo-colored */
  1364.     if (not_plane) get3_xy, xyzverts, xx, yy, zz, 1;
  1365.     else get3_xy, xyzverts, xx, yy;
  1366.     if (structof(vv)!=char)
  1367.       vv= split_bytscl(vv, 0, cmin=cmin, cmax=cmax);
  1368.   }
  1369.  
  1370.   /* accumulate nverts and values */
  1371.   item= numberof(nv)-1;
  1372.   nverts(list:list+item)= nv;
  1373.   values(list:list+item)= vv;
  1374.   list+= item+1;
  1375.  
  1376.   /* accumulate x, y, and z */
  1377.   item= numberof(xx)-1;
  1378.   x(vlist:vlist+item)= xx;
  1379.   y(vlist:vlist+item)= yy;
  1380.   if (not_plane) z(vlist:vlist+item)= zz;
  1381.   vlist+= item+1;
  1382. }
  1383.  
  1384. func _pl3tree_add(leaf, plane, tree)
  1385. {
  1386.   if (!is_void(_car(tree))) {
  1387.     /* tree has slicing plane, slice new leaf or plane and descend */
  1388.     back= _pl3tree_slice(_car(tree), leaf);
  1389.     if (back) {
  1390.       if (_car(tree,2)) _pl3tree_add, back, plane, _car(tree,2);
  1391.       else _car, tree, 2, _lst([], [], back, []);
  1392.     }
  1393.     if (leaf) {
  1394.       if (_car(tree,4)) _pl3tree_add, leaf, plane, _car(tree,4);
  1395.       else _car, tree, 4, _lst([], [], leaf, []);
  1396.     }
  1397.  
  1398.   } else if (!is_void(plane)) {
  1399.     /* tree is just a leaf, but this leaf has slicing plane */
  1400.     _car, tree, 1, plane;
  1401.     leaf= _car(tree, 3, leaf);  /* swap new leaf with original leaf */
  1402.     back= _pl3tree_slice(plane, leaf);   /* slice old leaf by plane */
  1403.     if (back) _car, tree, 2, _lst([], [], back, []);
  1404.     if (leaf) _car, tree, 4, _lst([], [], leaf, []);
  1405.  
  1406.   } else {
  1407.     /* tree is just a leaf and this leaf has no slicing plane */
  1408.     _cdr, leaf, 1, _car(tree, 3, leaf);
  1409.   }
  1410. }
  1411.  
  1412. func _pl3tree_slice(plane, &leaf)
  1413. {
  1414.   back= frnt= [];
  1415.   for ( ; leaf ; leaf=_cdr(leaf)) {
  1416.     ll= _car(leaf);
  1417.     /* each item in the leaf list is itself a list */
  1418.     nvf= nvb= _nxt(ll);
  1419.     xyzf= xyzb= _nxt(ll);
  1420.     valf= valb= _nxt(ll);
  1421.     slice2x, plane, nvf, xyzf, valf, nvb, xyzb, valb;
  1422.     if (numberof(nvf))
  1423.       frnt= _cat(_lst(_cat(nvf,xyzf,_lst(valf),ll)),frnt);
  1424.     if (numberof(nvb))
  1425.       back= _cat(_lst(_cat(nvb,xyzb,_lst(valb),ll)),back);
  1426.   }
  1427.   leaf= frnt;
  1428.   return back;
  1429. }
  1430.  
  1431. #if 0
  1432. func pl3tree_prt(void)
  1433. {
  1434.   tree= _cdr(_draw3_list, _draw3_n);
  1435.   if (is_void(tree) || _car(tree)!=pl3tree) {
  1436.     write, "<current 3D display not a pl3tree>";
  1437.   } else {
  1438.     tree= _car(tree,2);
  1439.     _pl3tree_prt,tree,0
  1440.   }
  1441. }
  1442.  
  1443. func _pl3tree_prt(tree,depth)
  1444. {
  1445.   if (is_void(tree)) return;
  1446.   indent= strpart("                               ",1:1+2*depth);
  1447.   indent= strpart(indent, 1:-1);
  1448.   write, indent+"+DEPTH= "+pr1(depth);
  1449.   if (_len(tree)!=4) write, indent+"***error - not a tree";
  1450.   write, indent+"plane= "+pr1(_car(tree));
  1451.   back= _car(tree,2);
  1452.   list= _car(tree,3);
  1453.   frnt= _car(tree,4);
  1454.   if (is_void(back)) write, indent+"back= []";
  1455.   else _pl3tree_prt, back, depth+1;
  1456.  
  1457.   while (list) {
  1458.     leaf= _nxt(list);
  1459.     write, indent+"leaf length= "+pr1(_len(leaf));
  1460.     write, indent+"npolys= "+pr1(numberof(_car(leaf)))+
  1461.       ", nverts= "+pr1(sum(_nxt(leaf)));
  1462.     write, indent+"nverts= "+pr1(numberof(_nxt(leaf))/3)+
  1463.       ", nvals= "+pr1(numberof(_nxt(leaf)));
  1464.   }
  1465.  
  1466.   if (is_void(frnt)) write, indent+"frnt= []";
  1467.   else _pl3tree_prt, frnt, depth+1;
  1468.   write, indent+"-DEPTH= "+pr1(depth);
  1469. }
  1470. #endif
  1471.  
  1472. /* ------------------------------------------------------------------------ */
  1473.